Número de SADs preditas não refutadas

Padrões gerais

Figura 1.1 número de SADs não refutadas ~ p * I(1-k) * MN. A linha azul é uma estimativa baseada em ‘loess’.

Para descrever estatisticamente o número de SADs preditas não refutadas iremos utilizar a distribuição binomial com função de ligação logito. Agruparemos os dados pelo sítio de amostragem. Esperamos que MNEE apresente melhor congruência com o observado para sítios em paisagens com proporção intermediária de cobertura vegetal e para cenários de extrema limitação de dispersão. Uma interpretação desta predição pode ser obtida se possibilitarmos que, para cada sítios de amostragem, a probabilidade de uma predição ser refutada seja uma função da variável de dispersão para cada modelo, ou seja, (dispersão * modelo neutro | sítio de amostragem). O modelo cheio considera a interação de terceira ordem entre as preditoras cobertura vegetal (p, contínua), dispersão (k, categorica; k_1, contínua; d_Lplot, contínua) e a classe do modelo neutro (MNEE e MNEI). Utilizamos uma abordagem baseada em seleção de modelos para determinar a relação entre as variáveis mais parcimoniosa.

GAMM binomial: número de predições não refutadas

A hipótese de trabalho é sobre uma relação não linear, uma escolha são generalized additive models onde há maior flexibilidade para estimar a tendência global dos dados (Wood 2017).

Comparação de Modelos Cheios

l_md <- vector("list",6)
names(l_md) <- c("k 1|SiteCode", "k MN|SiteCode",
                 "k_1 1|SiteCode","k_1 MN|SiteCode",
                 "d/L 1|SiteCode","d/L MN|SiteCode")
# muitos parametros para poucos dados "k_1 k_1 * MN|SiteCode","d_Lplot d_Lplot * MN|SiteCode")
# k
l_md[[1]] <- gam(cbind(n_nRef,100-n_nRef) ~ k * MN + s(p.z,k,by=MN,bs="fs") +
                   s(SiteCode,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
l_md[[2]] <- gam(cbind(n_nRef,100-n_nRef) ~ k * MN + s(p.z,k,by=MN,bs="fs") +
                   s(SiteCode,by=MN,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
# k_1
l_md[[3]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN + s(p.z,by=MN,bs="tp") + s(k_1.z,by=MN,bs="tp") + ti(p.z,k_1.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
l_md[[4]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN +  s(p.z,by=MN,bs="tp") + s(k_1.z,by=MN,bs="tp") + ti(p.z,k_1.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,by=MN,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
# d/L
l_md[[5]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN + s(p.z,by=MN,bs="tp") + s(d_Lplot.z,by=MN,bs="tp") + ti(p.z,d_Lplot.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
l_md[[6]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN +  s(p.z,by=MN,bs="tp") + s(d_Lplot.z,by=MN,bs="tp") + ti(p.z,d_Lplot.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,by=MN,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
AICctab(l_md,weights=TRUE)
##                 dAICc   df               weight
## k MN|SiteCode       0.0 567.902990229433 1     
## k_1 MN|SiteCode 10456.4 255.87941997766  <0.001
## d/L MN|SiteCode 14414.2 255.878481862603 <0.001
## k 1|SiteCode    29130.6 480.710534787585 <0.001
## k_1 1|SiteCode  38331.6 162.939249908097 <0.001
## d/L 1|SiteCode  43420.3 162.955609731983 <0.001
# save(l_md,file="~/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/Support Information/md_nRef.Rdata")
# load("~/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/Support Information/md_nRef.Rdata")

O único modelo plausível considera a variável k (fator) e estrutura aleatória MN|SiteCode.

Figura 1.2 appraise(gam( cbind(n_nRef,100-n_nRef) ~ k*MN + s(p.z,MN,by=k), binomial(logit) )

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## cbind(n_nRef, 100 - n_nRef) ~ k * MN + s(p.z, MN, by = k, bs = "fs") + 
##     s(SiteCode, by = MN, bs = "re")
## 
## Parametric coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   5.42767    0.91059   5.961 2.51e-09 ***
## k0.95         0.52243    1.12159   0.466 0.641365    
## k0.9          3.17125    1.43884   2.204 0.027523 *  
## k0.85        -1.87700    1.01948  -1.841 0.065602 .  
## k0.8         -2.98912    0.85346  -3.502 0.000461 ***
## k0.75        -2.92754    0.85370  -3.429 0.000605 ***
## k0.7         -2.95732    0.85436  -3.461 0.000537 ***
## k0.65        -2.94463    0.85909  -3.428 0.000609 ***
## k0.6         -2.74654    0.86452  -3.177 0.001488 ** 
## k0.55        -1.06282    1.40332  -0.757 0.448831    
## k0.5          7.28875    1.49804   4.866 1.14e-06 ***
## k0.45         7.49428    1.55653   4.815 1.47e-06 ***
## k0.4          2.79887    1.56415   1.789 0.073553 .  
## k0.35         5.34257    1.56681   3.410 0.000650 ***
## k0.3          6.95465    1.59726   4.354 1.34e-05 ***
## k0.25         5.88367    1.74370   3.374 0.000740 ***
## k0.2         10.84651    1.93458   5.607 2.06e-08 ***
## k0.15        10.33498    1.91456   5.398 6.74e-08 ***
## k0.1          9.21501    1.84974   4.982 6.30e-07 ***
## k0.05         1.96084    1.59389   1.230 0.218614    
## MNEI         -4.21846    1.12969  -3.734 0.000188 ***
## k0.95:MNEI   -1.14504    1.40522  -0.815 0.415159    
## k0.9:MNEI    -2.64414    1.76910  -1.495 0.135013    
## k0.85:MNEI    1.02849    1.28184   0.802 0.422351    
## k0.8:MNEI     2.23457    1.07546   2.078 0.037731 *  
## k0.75:MNEI    2.31624    1.07587   2.153 0.031326 *  
## k0.7:MNEI     2.53240    1.07688   2.352 0.018693 *  
## k0.65:MNEI    2.76399    1.08363   2.551 0.010752 *  
## k0.6:MNEI     2.67200    1.09096   2.449 0.014317 *  
## k0.55:MNEI    0.02772    1.69013   0.016 0.986914    
## k0.5:MNEI   -11.06077    1.88246  -5.876 4.21e-09 ***
## k0.45:MNEI  -10.17659    1.94165  -5.241 1.60e-07 ***
## k0.4:MNEI    -5.60540    1.94945  -2.875 0.004035 ** 
## k0.35:MNEI   -8.40553    1.95736  -4.294 1.75e-05 ***
## k0.3:MNEI   -11.12112    1.98664  -5.598 2.17e-08 ***
## k0.25:MNEI   -7.77766    2.10334  -3.698 0.000218 ***
## k0.2:MNEI   -14.38971    2.25911  -6.370 1.89e-10 ***
## k0.15:MNEI  -14.72112    2.23488  -6.587 4.49e-11 ***
## k0.1:MNEI   -15.59694    2.17908  -7.158 8.21e-13 ***
## k0.05:MNEI   -8.49448    1.99968  -4.248 2.16e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                        edf Ref.df    Chi.sq  p-value    
## s(p.z,MN):k0.99  1.186e+01     18  1316.247  < 2e-16 ***
## s(p.z,MN):k0.95  1.141e+01     18   709.588  < 2e-16 ***
## s(p.z,MN):k0.9   1.381e+01     18   444.347  < 2e-16 ***
## s(p.z,MN):k0.85  8.127e+00     18   208.742  < 2e-16 ***
## s(p.z,MN):k0.8   1.867e+00     18    63.259 9.46e-09 ***
## s(p.z,MN):k0.75  1.363e+00     18    11.048  0.00799 ** 
## s(p.z,MN):k0.7   4.283e-03     18     0.001  0.92382    
## s(p.z,MN):k0.65  2.511e+00     18    10.797  0.02563 *  
## s(p.z,MN):k0.6   3.271e+00     18    30.670 9.65e-05 ***
## s(p.z,MN):k0.55  1.156e+01     18   199.988  < 2e-16 ***
## s(p.z,MN):k0.5   1.652e+01     18  1595.015  < 2e-16 ***
## s(p.z,MN):k0.45  1.705e+01     18  2146.916  < 2e-16 ***
## s(p.z,MN):k0.4   1.743e+01     18  2889.126  < 2e-16 ***
## s(p.z,MN):k0.35  1.733e+01     18  2977.334  < 2e-16 ***
## s(p.z,MN):k0.3   1.754e+01     18  4704.837  < 2e-16 ***
## s(p.z,MN):k0.25  1.769e+01     18  6820.626  < 2e-16 ***
## s(p.z,MN):k0.2   1.775e+01     18  8094.560  < 2e-16 ***
## s(p.z,MN):k0.15  1.784e+01     18 14187.864  < 2e-16 ***
## s(p.z,MN):k0.1   1.786e+01     18 23904.904  < 2e-16 ***
## s(p.z,MN):k0.05  1.729e+01     18 42887.252  < 2e-16 ***
## s(SiteCode):MNEE 1.012e+02    102 34150.147  < 2e-16 ***
## s(SiteCode):MNEI 1.017e+02    102 32783.295  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.772   Deviance explained = 75.1%
## -REML =  41537  Scale est. = 1         n = 4120

Comparação com o observado

Figura 1.3 Observado e predito pelo modelo mais plausível. Eixo x = prop. de cobertura vegetal na paisagem; Eixo y: número de SADs não refutadas. Em branco o padrão geral estimado e em vermelho o erro padrão.

O padrão estimado parece ser muito sensível a valores extremo: a média de EI parece estar muito abaixo da mediana; essa tendência também parece ser observada para algum grau do vale que surge para p.z ~ -1 e p->0.05. Vejo duas opções: a) avaliar outliers e b) parâmetros da gam (parâmetro de penalidade, número de funções base; por padrão ambos são estimados utilizando REML, começar pelo parâmetro de penalidade)

Melhoria do modelo mais plausível: CONSTRUÇÃO

comparação de smoothers

comparação de parâmetro de penalização

diff_S = (S_obs - S_MN)/S_obs

Padrão Geral

Figura 2.1 diff_S = (S_obs - S_MN)/S_obs + (-1) * min(diff_S) + 0.01, inclinação_MNEE ~ 0 para todo k.

Parece que existem alguns outliers. De maneira geral MNEE apresenta boa congruência com os dados enquanto o comportamento de MNEI depende da dispersão: para k->100% o desvio aumento para ps extremos (bimodal/quadratica), apresentando boa congruência para proporção intermediárias de habitat, mas sempre subestimando S_obs. Com o aumento da capacidade de dispersão as estimavas de S se tornam mais próximas ao observado para p>0.5, superestimando o observado para valores elevados de p e k>0.5 o que leva a reduzir a porção de p em que MNEI faz uma boa aproximação.

GAMM: diff_S

Para modelar o padrão não linear dos dados vou utilizar GAMM: normal, Gamma(log) e Gamma(inverse) com duas estruturas aleatórias (1|SiteCode e MN|SiteCode).

Modelo Cheio e Gráficos Diagnosticos

##                           dAICc  df               weight
## normal id MN|SiteCode        0.0 204.716290217321 1     
## normal id 1|SiteCode      3500.4 201.919668132722 <0.001
## Gamma inverse MN|SiteCode 7545.2 183.096851367051 <0.001
## Gamma log MN|SiteCode     7899.2 183.888877817823 <0.001
## Gamma inverse 1|SiteCode  8599.7 172.829381805835 <0.001
## Gamma log 1|SiteCode      8730.6 169.279594585853 <0.001

Figura 2.2 Graficos diagnosticos do modelo mais plausível para diff_S.

Parece que o pressuposto de distribuição normal não é adequado, assim como parece haver alguma tendência nos dados que não está sendo detectada. Especulo que existam outliers que estão inflenciando muito a estimativa dos dados.

Observado e Predito

Figura 2.3 Diferença na riqueza observada e estimada pela proporção de cobertura vegetal.

Os modelos ajustados para MNEI estão sistematicamente com problemas de errar a média (figura 1.1 e 1.4; figura 2.1 e 2.3). Especulo que talvez alguns pontos estão influenciando muito as estimativas centrais ou estou usando configuração errada nas funções. Além disso para MNEI percebo que a tendência geral não está descrevendo o padrão dos dados, em especial para k->1.

U - taxa de especiação necessária para obter a riqueza observada no equilíbrio

Padrões Gerais

Figura 3.1 Padrões gerais de U_med: ~ k (% de propágulos até a vizinhança imediata da planta progenitora); ~ Stotal (riqueza observada na área amostral). E o padrão empírico encontrado nos dados Stotal ~ p

GAMM

Como os modelos com estrutura aleatória mais complexa levam muito tempo para estimar, irei começar com a estrutura aleatória mais simples para avaliar qual a melhor variável de dispersão e se for necessário refazer a seleção da estrutura aleatória.

Tabela de Seleção do Modelo Mais plausível:

##                    dAICc  df               weight
## d log gamma           0.0 121.264568083601 1     
## 1-k log gamma       174.1 125.474107334272 <0.001
## k.f log gamma       227.4 154.787008201268 <0.001
## d log normal        394.5 124.471358521882 <0.001
## d inverse gamma     766.2 127.68541117751  <0.001
## 1-k log normal      855.7 127.94024990325  <0.001
## k.f log normal      889.4 170.667471702588 <0.001
## 1-k inverse gamma   983.2 124.456458934862 <0.001
## k.f inverse gamma  1002.8 140.80259373167  <0.001
## d inverse normal   1136.4 121.957626137592 <0.001
## k.f inverse normal 1345.5 165.620770162204 <0.001
## 1-k inverse normal 1354.2 122.442084920132 <0.001
## 1-k id normal      2163.0 125.186614559034 <0.001
## d id normal        2189.1 123.301338427342 <0.001
## k.f id normal      2242.8 160.04673943848  <0.001

Segue gráficos diagnosticos do modelo mais plausível

Figura 3.2 Gráficos Diagnóstico do modelo mais plausível

O modelo apresenta boa congruência com os dados, porém parece existir outliers no entando.

Figura 3.3 Métricas de influência dos dados (broom::augment())

Exclui os dois sítios com valores de grande influência do conjunto de dados:

Figura 3.4 Grafico diagnostico do modelo mais plausivel sem os outliers.

Figura 4.5 Efeitos Parciais do modelo mais plausível

Parâmetros dos Modelos Neutros

d / L_plot

padrões gerais

Figura 5.1.1 d/Lplot ~ k

Descrição Estatística

## [1] "exp(coef):"
## (Intercept)       k0.95        k0.9       k0.85        k0.8       k0.75 
##   0.6715478   1.4153327   1.7360273   2.0081002   2.2669913   2.5268653 
##        k0.7       k0.65        k0.6       k0.55        k0.5       k0.45 
##   2.7936350   3.0726860   3.3755022   3.7066527   4.0756037   4.4953275 
##        k0.4       k0.35        k0.3       k0.25        k0.2       k0.15 
##   4.9791588   5.5515868   6.2545191   7.1436593   8.3360091  10.0630874 
##        k0.1       k0.05 
##  12.9390064  19.3670608

Figura 5.1.2 Gráficos Diagnósticos do modelo gam(d ~ k + s(DA,by=k), family=Gamma(link=“log”)). Exponenciação dos coeficientes parametricos estimados pelo modelo.

m’

padrões gerais

Figura 5.2.1 m’ ~ p (~k). A variável m’ é uma função de p, d e I(1/J)

Descrição Estatística

##        dAICc  df               weight
## model2    0.0 158.703955395811 1     
## model3  488.0 43               <0.001
## model1 7043.7 113.502207953068 <0.001

theta

padrões gerais

Descrição Estatística

–>

–> –> –> –>

–>

–> –> –> –>

–> –> –> –>

–>

–>

–>

–>

–>

–> –> –> –> –>

–> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–>

–> –> –> –>

–>

–> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–> –> –> –>

–>

–>

–>

–> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>